clc
clear all

n = 500;
% m = 30; % order of cyclic group
p = 0.3;   % observation probability
q = 0.3; % probability of not being corrupted by (multiplicative) outlier noise, q=1 means no outlier noise

Number_of_m = 11;   % number of points
Array_m = 6 + 4 * (1:Number_of_m)';   


Number_of_Sim = 30;   % number of repeats

% error 
MCR_Spectral = zeros(Number_of_Sim, Number_of_m);
MCR_Spectral_K = zeros(Number_of_Sim, Number_of_m);
MCR_PPM = zeros(Number_of_Sim, Number_of_m);
MCR_GPM = zeros(Number_of_Sim, Number_of_m);

% running time
time_Spectral = zeros(Number_of_Sim, Number_of_m);
time_Spectral_K = zeros(Number_of_Sim, Number_of_m);
time_PPM = zeros(Number_of_Sim, Number_of_m);
time_GPM = zeros(Number_of_Sim, Number_of_m);

for index_m = 1:Number_of_m
    m = Array_m(index_m);
    K = 10;


    
    %% Generating instance...
    % Set parameters
    param.n             = n;     % number of variables
    param.m             = m;      % number of states
    param.p_obs         = p;       % observation ratio
    param.maxIter       = 50;      % maximum number of projected power iterations
    param.mu0           = 10;      % scale factor
    param.pi0           = q;     % non-corruption rate
    param.round_flag    = 'no';    % indicate whether mu0 is infinity
    param.max_poweriter = 20;      % maximum number of power iterations
    param.eps_power     = 0.02;    % tolerance level of power methods

    % Noise density for the random corruption model
    P0          = (1-param.pi0)* ones( param.m, 1) / param.m;
    P0(1)       = P0(1) + param.pi0;
    param.P0    = P0;
    % display(sprintf('n=%d, m=%d, corruption_rate=%f', param.n, param.m, 1-param.pi0));

    % Generate a synthetic example
    parfor index_sim = 1:Number_of_Sim

        fprintf('index_n: %d of %d, index_sim: %d of %d\n', index_m, Number_of_m, index_sim, Number_of_Sim);

        [Problem, C, S_true]     = ProblemGenerator(param);


        %% Projected power method
        tic;
        L           = Problem.L;        % Input log-likelihood matrix
        Lambda      = eigs( L, 2);      % Compute 2nd largest eigenvalue of L
        L_opt       = @(I) L  * I;      % Matrix-vector multiplication operator
        proj        = @(z) ProjectOperator(z, param.n, param.m, param.mu0/Lambda(2)); % Projection operator
        [z, relErrors]  = noncvxAlign(param, Problem.X_gt, L_opt, proj);  % z is the output and relErrors is the L2 error
        t_PPM = toc;
        S_PPM = z2S(z, m);

        %% Spectral Estimator
        [t_Spectral, S_Spectral] = Spectral_Cyclic_old(C, m);

        [t_Spectral_K, S_Spectral_K] = Spectral_Cyclic(C, m, K);

        %% GPM
        [t_GPM, S_GPM] = GPM_Cyclic( C, S_Spectral_K, m );


        %% Update error
        MCR_Spectral(index_sim, index_m) = MCR(S_true, S_Spectral, m);
        MCR_Spectral_K(index_sim, index_m) = MCR(S_true, S_Spectral_K, m);
        MCR_GPM(index_sim, index_m) = MCR(S_true, S_GPM, m);
        MCR_PPM(index_sim, index_m) = MCR(S_true, S_PPM, m);

        %% Update time
        time_Spectral(index_sim, index_m) = t_Spectral;
        time_Spectral_K(index_sim, index_m) = t_Spectral_K;
        time_GPM(index_sim, index_m) = t_GPM + t_Spectral;
        time_PPM(index_sim, index_m) =  t_PPM;
    end
end

% 
fig = figure;
set(gcf,'units','points','position',[10,10,800,400])
% title('d=5, observe p_1 = 0.1, noise p_2=0.7')

subplot(1,2,1);
set(gca,'fontsize',15)
plot(Array_m, 1 - sum(MCR_PPM)/Number_of_Sim,'-p','LineWidth',2,'MarkerSize',10);
hold on
plot(Array_m, 1 - sum(MCR_GPM)/Number_of_Sim,'-o','LineWidth',2,'MarkerSize',10);
plot(Array_m, 1 - sum(MCR_Spectral_K)/Number_of_Sim,'-p','LineWidth',2,'MarkerSize',10);
plot(Array_m, 1 - sum(MCR_Spectral)/Number_of_Sim,'-*','LineWidth',2,'MarkerSize',10);
hold off
legend('PPM','GPM','Spectral-K', 'Spectral','Location','southwest')
xlabel('$m$','Interpreter','latex')
ylabel('\textbf{Recovery Rate}','Interpreter','latex')


subplot(1,2,2);
set(gca,'fontsize',15)
plot(Array_m,sum(time_PPM)/(Number_of_Sim),'-p','LineWidth',2,'MarkerSize',10);
hold on
plot(Array_m,sum(time_GPM)/(Number_of_Sim),'-o','LineWidth',2,'MarkerSize',10);
plot(Array_m,sum(time_Spectral_K)/(Number_of_Sim),'-p','LineWidth',2,'MarkerSize',10);
plot(Array_m,sum(time_Spectral)/(Number_of_Sim),'-p','LineWidth',2,'MarkerSize',10);
hold off
legend('PPM','GPM','Spectral-K','Spectral', 'Location','northwest')
xlabel('$m$','Interpreter','latex')
ylabel('\textbf{Time (s)}','Interpreter','latex')

% sgtitle(['$n$ = ', num2str(n), ', $p$ = ', num2str(p), ', $q$ = ', num2str(q)],'Interpreter','latex');

% Save
folder_name = ['n=', num2str(n),'_m=',num2str(Array_m(1)),'-',num2str(Array_m(Number_of_m)), '_p=', strrep(num2str(p),'.',''), '_q=', strrep(num2str(q),'.','')];
mkdir('C:\Users\Administrator\Desktop\Cyclic\Result',folder_name);
save(['Result\',num2str(folder_name),'\MCR_Spec.mat'], 'MCR_Spectral');
save(['Result\',num2str(folder_name),'\MCR_Spec_K.mat'], 'MCR_Spectral_K');
save(['Result\',num2str(folder_name),'\MCR_GPM.mat'], 'MCR_GPM');
save(['Result\',num2str(folder_name),'\MCR_PPM.mat'], 'MCR_PPM');
save(['Result\',num2str(folder_name),'\time_Spec.mat'], 'time_Spectral');
save(['Result\',num2str(folder_name),'\time_Spec_K.mat'], 'time_Spectral_K');
save(['Result\',num2str(folder_name),'\time_GPM.mat'], 'time_GPM');
save(['Result\',num2str(folder_name),'\time_PPM.mat'], 'time_PPM');
savefig(['Result\',num2str(folder_name),'\',num2str(folder_name),'.fig']);

clc
clear all

n = 500;
% m = 30; % order of cyclic group
p = 0.3;   % observation probability
q = 0.7; % probability of not being corrupted by (multiplicative) outlier noise, q=1 means no outlier noise

Number_of_m = 11;   % number of points
Array_m = 6 + 4 * (1:Number_of_m)';   


Number_of_Sim = 30;   % number of repeats

% error 
MCR_Spectral = zeros(Number_of_Sim, Number_of_m);
MCR_Spectral_K = zeros(Number_of_Sim, Number_of_m);
MCR_PPM = zeros(Number_of_Sim, Number_of_m);
MCR_GPM = zeros(Number_of_Sim, Number_of_m);

% running time
time_Spectral = zeros(Number_of_Sim, Number_of_m);
time_Spectral_K = zeros(Number_of_Sim, Number_of_m);
time_PPM = zeros(Number_of_Sim, Number_of_m);
time_GPM = zeros(Number_of_Sim, Number_of_m);

for index_m = 1:Number_of_m
    m = Array_m(index_m);
    K = 10;


    
    %% Generating instance...
    % Set parameters
    param.n             = n;     % number of variables
    param.m             = m;      % number of states
    param.p_obs         = p;       % observation ratio
    param.maxIter       = 50;      % maximum number of projected power iterations
    param.mu0           = 10;      % scale factor
    param.pi0           = q;     % non-corruption rate
    param.round_flag    = 'no';    % indicate whether mu0 is infinity
    param.max_poweriter = 20;      % maximum number of power iterations
    param.eps_power     = 0.02;    % tolerance level of power methods

    % Noise density for the random corruption model
    P0          = (1-param.pi0)* ones( param.m, 1) / param.m;
    P0(1)       = P0(1) + param.pi0;
    param.P0    = P0;
    % display(sprintf('n=%d, m=%d, corruption_rate=%f', param.n, param.m, 1-param.pi0));

    % Generate a synthetic example
    parfor index_sim = 1:Number_of_Sim

        fprintf('index_n: %d of %d, index_sim: %d of %d\n', index_m, Number_of_m, index_sim, Number_of_Sim);

        [Problem, C, S_true]     = ProblemGenerator(param);


        %% Projected power method
        tic;
        L           = Problem.L;        % Input log-likelihood matrix
        Lambda      = eigs( L, 2);      % Compute 2nd largest eigenvalue of L
        L_opt       = @(I) L  * I;      % Matrix-vector multiplication operator
        proj        = @(z) ProjectOperator(z, param.n, param.m, param.mu0/Lambda(2)); % Projection operator
        [z, relErrors]  = noncvxAlign(param, Problem.X_gt, L_opt, proj);  % z is the output and relErrors is the L2 error
        t_PPM = toc;
        S_PPM = z2S(z, m);

        %% Spectral Estimator
        [t_Spectral, S_Spectral] = Spectral_Cyclic_old(C, m);

        [t_Spectral_K, S_Spectral_K] = Spectral_Cyclic(C, m, K);

        %% GPM
        [t_GPM, S_GPM] = GPM_Cyclic( C, S_Spectral_K, m );


        %% Update error
        MCR_Spectral(index_sim, index_m) = MCR(S_true, S_Spectral, m);
        MCR_Spectral_K(index_sim, index_m) = MCR(S_true, S_Spectral_K, m);
        MCR_GPM(index_sim, index_m) = MCR(S_true, S_GPM, m);
        MCR_PPM(index_sim, index_m) = MCR(S_true, S_PPM, m);

        %% Update time
        time_Spectral(index_sim, index_m) = t_Spectral;
        time_Spectral_K(index_sim, index_m) = t_Spectral_K;
        time_GPM(index_sim, index_m) = t_GPM + t_Spectral;
        time_PPM(index_sim, index_m) =  t_PPM;
    end
end

% 
fig = figure;
set(gcf,'units','points','position',[10,10,800,400])
% title('d=5, observe p_1 = 0.1, noise p_2=0.7')

subplot(1,2,1);
set(gca,'fontsize',15)
plot(Array_m, 1 - sum(MCR_PPM)/Number_of_Sim,'-p','LineWidth',2,'MarkerSize',10);
hold on
plot(Array_m, 1 - sum(MCR_GPM)/Number_of_Sim,'-o','LineWidth',2,'MarkerSize',10);
plot(Array_m, 1 - sum(MCR_Spectral_K)/Number_of_Sim,'-p','LineWidth',2,'MarkerSize',10);
plot(Array_m, 1 - sum(MCR_Spectral)/Number_of_Sim,'-*','LineWidth',2,'MarkerSize',10);
hold off
legend('PPM','GPM','Spectral-K', 'Spectral','Location','southwest')
xlabel('$m$','Interpreter','latex')
ylabel('\textbf{Recovery Rate}','Interpreter','latex')


subplot(1,2,2);
set(gca,'fontsize',15)
plot(Array_m,sum(time_PPM)/(Number_of_Sim),'-p','LineWidth',2,'MarkerSize',10);
hold on
plot(Array_m,sum(time_GPM)/(Number_of_Sim),'-o','LineWidth',2,'MarkerSize',10);
plot(Array_m,sum(time_Spectral_K)/(Number_of_Sim),'-p','LineWidth',2,'MarkerSize',10);
plot(Array_m,sum(time_Spectral)/(Number_of_Sim),'-p','LineWidth',2,'MarkerSize',10);
hold off
legend('PPM','GPM','Spectral-K','Spectral', 'Location','northwest')
xlabel('$m$','Interpreter','latex')
ylabel('\textbf{Time (s)}','Interpreter','latex')

% sgtitle(['$n$ = ', num2str(n), ', $p$ = ', num2str(p), ', $q$ = ', num2str(q)],'Interpreter','latex');

% Save
folder_name = ['n=', num2str(n),'_m=',num2str(Array_m(1)),'-',num2str(Array_m(Number_of_m)), '_p=', strrep(num2str(p),'.',''), '_q=', strrep(num2str(q),'.','')];
mkdir('C:\Users\Administrator\Desktop\Cyclic\Result',folder_name);
save(['Result\',num2str(folder_name),'\MCR_Spec.mat'], 'MCR_Spectral');
save(['Result\',num2str(folder_name),'\MCR_Spec_K.mat'], 'MCR_Spectral_K');
save(['Result\',num2str(folder_name),'\MCR_GPM.mat'], 'MCR_GPM');
save(['Result\',num2str(folder_name),'\MCR_PPM.mat'], 'MCR_PPM');
save(['Result\',num2str(folder_name),'\time_Spec.mat'], 'time_Spectral');
save(['Result\',num2str(folder_name),'\time_Spec_K.mat'], 'time_Spectral_K');
save(['Result\',num2str(folder_name),'\time_GPM.mat'], 'time_GPM');
save(['Result\',num2str(folder_name),'\time_PPM.mat'], 'time_PPM');
savefig(['Result\',num2str(folder_name),'\',num2str(folder_name),'.fig']);


clc
clear all

n = 500;
% m = 30; % order of cyclic group
p = 0.7;   % observation probability
q = 0.3; % probability of not being corrupted by (multiplicative) outlier noise, q=1 means no outlier noise

Number_of_m = 11;   % number of points
Array_m = 6 + 4 * (1:Number_of_m)';   


Number_of_Sim = 30;   % number of repeats

% error 
MCR_Spectral = zeros(Number_of_Sim, Number_of_m);
MCR_Spectral_K = zeros(Number_of_Sim, Number_of_m);
MCR_PPM = zeros(Number_of_Sim, Number_of_m);
MCR_GPM = zeros(Number_of_Sim, Number_of_m);

% running time
time_Spectral = zeros(Number_of_Sim, Number_of_m);
time_Spectral_K = zeros(Number_of_Sim, Number_of_m);
time_PPM = zeros(Number_of_Sim, Number_of_m);
time_GPM = zeros(Number_of_Sim, Number_of_m);

for index_m = 1:Number_of_m
    m = Array_m(index_m);
    K = 10;


    
    %% Generating instance...
    % Set parameters
    param.n             = n;     % number of variables
    param.m             = m;      % number of states
    param.p_obs         = p;       % observation ratio
    param.maxIter       = 50;      % maximum number of projected power iterations
    param.mu0           = 10;      % scale factor
    param.pi0           = q;     % non-corruption rate
    param.round_flag    = 'no';    % indicate whether mu0 is infinity
    param.max_poweriter = 20;      % maximum number of power iterations
    param.eps_power     = 0.02;    % tolerance level of power methods

    % Noise density for the random corruption model
    P0          = (1-param.pi0)* ones( param.m, 1) / param.m;
    P0(1)       = P0(1) + param.pi0;
    param.P0    = P0;
    % display(sprintf('n=%d, m=%d, corruption_rate=%f', param.n, param.m, 1-param.pi0));

    % Generate a synthetic example
    parfor index_sim = 1:Number_of_Sim

        fprintf('index_n: %d of %d, index_sim: %d of %d\n', index_m, Number_of_m, index_sim, Number_of_Sim);

        [Problem, C, S_true]     = ProblemGenerator(param);


        %% Projected power method
        tic;
        L           = Problem.L;        % Input log-likelihood matrix
        Lambda      = eigs( L, 2);      % Compute 2nd largest eigenvalue of L
        L_opt       = @(I) L  * I;      % Matrix-vector multiplication operator
        proj        = @(z) ProjectOperator(z, param.n, param.m, param.mu0/Lambda(2)); % Projection operator
        [z, relErrors]  = noncvxAlign(param, Problem.X_gt, L_opt, proj);  % z is the output and relErrors is the L2 error
        t_PPM = toc;
        S_PPM = z2S(z, m);

        %% Spectral Estimator
        [t_Spectral, S_Spectral] = Spectral_Cyclic_old(C, m);

        [t_Spectral_K, S_Spectral_K] = Spectral_Cyclic(C, m, K);

        %% GPM
        [t_GPM, S_GPM] = GPM_Cyclic( C, S_Spectral_K, m );


        %% Update error
        MCR_Spectral(index_sim, index_m) = MCR(S_true, S_Spectral, m);
        MCR_Spectral_K(index_sim, index_m) = MCR(S_true, S_Spectral_K, m);
        MCR_GPM(index_sim, index_m) = MCR(S_true, S_GPM, m);
        MCR_PPM(index_sim, index_m) = MCR(S_true, S_PPM, m);

        %% Update time
        time_Spectral(index_sim, index_m) = t_Spectral;
        time_Spectral_K(index_sim, index_m) = t_Spectral_K;
        time_GPM(index_sim, index_m) = t_GPM + t_Spectral;
        time_PPM(index_sim, index_m) =  t_PPM;
    end
end

% 
fig = figure;
set(gcf,'units','points','position',[10,10,800,400])
% title('d=5, observe p_1 = 0.1, noise p_2=0.7')

subplot(1,2,1);
set(gca,'fontsize',15)
plot(Array_m, 1 - sum(MCR_PPM)/Number_of_Sim,'-p','LineWidth',2,'MarkerSize',10);
hold on
plot(Array_m, 1 - sum(MCR_GPM)/Number_of_Sim,'-o','LineWidth',2,'MarkerSize',10);
plot(Array_m, 1 - sum(MCR_Spectral_K)/Number_of_Sim,'-p','LineWidth',2,'MarkerSize',10);
plot(Array_m, 1 - sum(MCR_Spectral)/Number_of_Sim,'-*','LineWidth',2,'MarkerSize',10);
hold off
legend('PPM','GPM','Spectral-K', 'Spectral','Location','southwest')
xlabel('$m$','Interpreter','latex')
ylabel('\textbf{Recovery Rate}','Interpreter','latex')


subplot(1,2,2);
set(gca,'fontsize',15)
plot(Array_m,sum(time_PPM)/(Number_of_Sim),'-p','LineWidth',2,'MarkerSize',10);
hold on
plot(Array_m,sum(time_GPM)/(Number_of_Sim),'-o','LineWidth',2,'MarkerSize',10);
plot(Array_m,sum(time_Spectral_K)/(Number_of_Sim),'-p','LineWidth',2,'MarkerSize',10);
plot(Array_m,sum(time_Spectral)/(Number_of_Sim),'-p','LineWidth',2,'MarkerSize',10);
hold off
legend('PPM','GPM','Spectral-K','Spectral', 'Location','northwest')
xlabel('$m$','Interpreter','latex')
ylabel('\textbf{Time (s)}','Interpreter','latex')

% sgtitle(['$n$ = ', num2str(n), ', $p$ = ', num2str(p), ', $q$ = ', num2str(q)],'Interpreter','latex');

% Save
folder_name = ['n=', num2str(n),'_m=',num2str(Array_m(1)),'-',num2str(Array_m(Number_of_m)), '_p=', strrep(num2str(p),'.',''), '_q=', strrep(num2str(q),'.','')];
mkdir('C:\Users\Administrator\Desktop\Cyclic\Result',folder_name);
save(['Result\',num2str(folder_name),'\MCR_Spec.mat'], 'MCR_Spectral');
save(['Result\',num2str(folder_name),'\MCR_Spec_K.mat'], 'MCR_Spectral_K');
save(['Result\',num2str(folder_name),'\MCR_GPM.mat'], 'MCR_GPM');
save(['Result\',num2str(folder_name),'\MCR_PPM.mat'], 'MCR_PPM');
save(['Result\',num2str(folder_name),'\time_Spec.mat'], 'time_Spectral');
save(['Result\',num2str(folder_name),'\time_Spec_K.mat'], 'time_Spectral_K');
save(['Result\',num2str(folder_name),'\time_GPM.mat'], 'time_GPM');
save(['Result\',num2str(folder_name),'\time_PPM.mat'], 'time_PPM');
savefig(['Result\',num2str(folder_name),'\',num2str(folder_name),'.fig']);



clc
clear all

n = 500;
% m = 30; % order of cyclic group
p = 0.7;   % observation probability
q = 0.7; % probability of not being corrupted by (multiplicative) outlier noise, q=1 means no outlier noise

Number_of_m = 11;   % number of points
Array_m = 6 + 4 * (1:Number_of_m)';   


Number_of_Sim = 30;   % number of repeats

% error 
MCR_Spectral = zeros(Number_of_Sim, Number_of_m);
MCR_Spectral_K = zeros(Number_of_Sim, Number_of_m);
MCR_PPM = zeros(Number_of_Sim, Number_of_m);
MCR_GPM = zeros(Number_of_Sim, Number_of_m);

% running time
time_Spectral = zeros(Number_of_Sim, Number_of_m);
time_Spectral_K = zeros(Number_of_Sim, Number_of_m);
time_PPM = zeros(Number_of_Sim, Number_of_m);
time_GPM = zeros(Number_of_Sim, Number_of_m);

for index_m = 1:Number_of_m
    m = Array_m(index_m);
    K = 10;


    
    %% Generating instance...
    % Set parameters
    param.n             = n;     % number of variables
    param.m             = m;      % number of states
    param.p_obs         = p;       % observation ratio
    param.maxIter       = 50;      % maximum number of projected power iterations
    param.mu0           = 10;      % scale factor
    param.pi0           = q;     % non-corruption rate
    param.round_flag    = 'no';    % indicate whether mu0 is infinity
    param.max_poweriter = 20;      % maximum number of power iterations
    param.eps_power     = 0.02;    % tolerance level of power methods

    % Noise density for the random corruption model
    P0          = (1-param.pi0)* ones( param.m, 1) / param.m;
    P0(1)       = P0(1) + param.pi0;
    param.P0    = P0;
    % display(sprintf('n=%d, m=%d, corruption_rate=%f', param.n, param.m, 1-param.pi0));

    % Generate a synthetic example
    parfor index_sim = 1:Number_of_Sim

        fprintf('index_n: %d of %d, index_sim: %d of %d\n', index_m, Number_of_m, index_sim, Number_of_Sim);

        [Problem, C, S_true]     = ProblemGenerator(param);


        %% Projected power method
        tic;
        L           = Problem.L;        % Input log-likelihood matrix
        Lambda      = eigs( L, 2);      % Compute 2nd largest eigenvalue of L
        L_opt       = @(I) L  * I;      % Matrix-vector multiplication operator
        proj        = @(z) ProjectOperator(z, param.n, param.m, param.mu0/Lambda(2)); % Projection operator
        [z, relErrors]  = noncvxAlign(param, Problem.X_gt, L_opt, proj);  % z is the output and relErrors is the L2 error
        t_PPM = toc;
        S_PPM = z2S(z, m);

        %% Spectral Estimator
        [t_Spectral, S_Spectral] = Spectral_Cyclic_old(C, m);

        [t_Spectral_K, S_Spectral_K] = Spectral_Cyclic(C, m, K);

        %% GPM
        [t_GPM, S_GPM] = GPM_Cyclic( C, S_Spectral_K, m );


        %% Update error
        MCR_Spectral(index_sim, index_m) = MCR(S_true, S_Spectral, m);
        MCR_Spectral_K(index_sim, index_m) = MCR(S_true, S_Spectral_K, m);
        MCR_GPM(index_sim, index_m) = MCR(S_true, S_GPM, m);
        MCR_PPM(index_sim, index_m) = MCR(S_true, S_PPM, m);

        %% Update time
        time_Spectral(index_sim, index_m) = t_Spectral;
        time_Spectral_K(index_sim, index_m) = t_Spectral_K;
        time_GPM(index_sim, index_m) = t_GPM + t_Spectral;
        time_PPM(index_sim, index_m) =  t_PPM;
    end
end

% 
fig = figure;
set(gcf,'units','points','position',[10,10,800,400])
% title('d=5, observe p_1 = 0.1, noise p_2=0.7')

subplot(1,2,1);
set(gca,'fontsize',15)
plot(Array_m, 1 - sum(MCR_PPM)/Number_of_Sim,'-p','LineWidth',2,'MarkerSize',10);
hold on
plot(Array_m, 1 - sum(MCR_GPM)/Number_of_Sim,'-o','LineWidth',2,'MarkerSize',10);
plot(Array_m, 1 - sum(MCR_Spectral_K)/Number_of_Sim,'-p','LineWidth',2,'MarkerSize',10);
plot(Array_m, 1 - sum(MCR_Spectral)/Number_of_Sim,'-*','LineWidth',2,'MarkerSize',10);
hold off
legend('PPM','GPM','Spectral-K', 'Spectral','Location','southwest')
xlabel('$m$','Interpreter','latex')
ylabel('\textbf{Recovery Rate}','Interpreter','latex')


subplot(1,2,2);
set(gca,'fontsize',15)
plot(Array_m,sum(time_PPM)/(Number_of_Sim),'-p','LineWidth',2,'MarkerSize',10);
hold on
plot(Array_m,sum(time_GPM)/(Number_of_Sim),'-o','LineWidth',2,'MarkerSize',10);
plot(Array_m,sum(time_Spectral_K)/(Number_of_Sim),'-p','LineWidth',2,'MarkerSize',10);
plot(Array_m,sum(time_Spectral)/(Number_of_Sim),'-p','LineWidth',2,'MarkerSize',10);
hold off
legend('PPM','GPM','Spectral-K','Spectral', 'Location','northwest')
xlabel('$m$','Interpreter','latex')
ylabel('\textbf{Time (s)}','Interpreter','latex')

% sgtitle(['$n$ = ', num2str(n), ', $p$ = ', num2str(p), ', $q$ = ', num2str(q)],'Interpreter','latex');

% Save
folder_name = ['n=', num2str(n),'_m=',num2str(Array_m(1)),'-',num2str(Array_m(Number_of_m)), '_p=', strrep(num2str(p),'.',''), '_q=', strrep(num2str(q),'.','')];
mkdir('C:\Users\Administrator\Desktop\Cyclic\Result',folder_name);
save(['Result\',num2str(folder_name),'\MCR_Spec.mat'], 'MCR_Spectral');
save(['Result\',num2str(folder_name),'\MCR_Spec_K.mat'], 'MCR_Spectral_K');
save(['Result\',num2str(folder_name),'\MCR_GPM.mat'], 'MCR_GPM');
save(['Result\',num2str(folder_name),'\MCR_PPM.mat'], 'MCR_PPM');
save(['Result\',num2str(folder_name),'\time_Spec.mat'], 'time_Spectral');
save(['Result\',num2str(folder_name),'\time_Spec_K.mat'], 'time_Spectral_K');
save(['Result\',num2str(folder_name),'\time_GPM.mat'], 'time_GPM');
save(['Result\',num2str(folder_name),'\time_PPM.mat'], 'time_PPM');
savefig(['Result\',num2str(folder_name),'\',num2str(folder_name),'.fig']);

